/*
 * Wells.cpp
 *
 *  Created on: 21/09/2012
 *      Author: rogsoares
 */

#include "SimulatorParameters.h"

namespace PRS{

// get node ids associated to wells
void SimulatorParameters::getWells(TMesh* pTMesh_data){
	/*
	 * Loop over all nodes and boundary elements searching for well flags
	 * Store them on a temporary set<VertexData*>
	 */
	int dim = pTMesh_data->getMeshDim();
	VertexData* vdata;
	BdryElementData *belem;
	typedef std::set<VertexData*> setWells;
	std::map<int,setWells> mapPWells;
	std::map<int,setWells> mapIWells;
	int ndom = pTMesh_data->getNumDomains();
	for(int i=0; i<ndom; i++){
		int nrows, ncols;
		/*
		 * Loop over all nodes
		 */
		pTMesh_data->vertices_list[i].getSize(nrows,ncols);
		for(int j=0; j<nrows; j++){
			vdata = pTMesh_data->vertices_list[i].getValue(j,0);
			int flag = vdata->flag;
			if ( this->isProductionWell(flag) ){
				mapPWells[flag].insert(vdata);
			}
			if ( this->isInjectionWell(flag) ){
				mapIWells[flag].insert(vdata);
			}
		}
		/*
		 * loop over boundary elements
		 */
		pTMesh_data->bdry_elem_list[i].getSize(nrows,ncols);
		for(int j=0; j<nrows; j++){
			belem = pTMesh_data->bdry_elem_list[i].getValue(j,0);
			VertexData *vecdata[3] = {belem->vertex_0,belem->vertex_1,belem->vertex_2};
			for (int k=0;k<dim;k++){
				int flag = vecdata[k]->flag;
				if ( this->isProductionWell(flag) ){
					mapPWells[flag].insert(vecdata[k]);
				}
				if ( this->isInjectionWell(flag) ){
					mapIWells[flag].insert(vecdata[k]);
				}
			}
		}
	}

	/*
	 * Once all nodes flagged as well has been detected, let's organize them as:
	 * Injection wells:
	 * flag_1 : list of nodes
	 * flag_2 : list of nodes
	 * ...
	 * flag_n : list of nodes
	 *
	 * Production wells:
	 * idem Injection wells
	 */

	int i=0;
	VertexWell *vwell;
	int numPW = mapPWells.size();												// total number of production wells
	pTMesh_data->setNumProductionWells(numPW);
	pTMesh_data->productionWells_list = new Matrix<VertexWell*>[numPW];			// allocate memory for a set of node well for each flag
	// loop over all nodes flagged as injection wells
	std::map<int,setWells>::iterator mapIter = mapPWells.begin();
	for(;mapIter!=mapPWells.end();mapIter++){
		int flag = mapIter->first;
		pTMesh_data->setProductionWellFlag(i,flag);
		well_info = MWells[flag];
		double Qt = well_info.flowRate;
		double vt = .0;
		setWells setPW = mapIter->second;
		int num =  setPW.size();												// number of node for well flag
		pTMesh_data->productionWells_list[i].allocateMemory(num,1);
		// caculate total volume for a specific well. Flow rate for all nodes of the specific flag must be weighted: Qi = Qt(Volume_node/Total_Volume)
		setWells::iterator setIter = setPW.begin();
		for(;setIter!=setPW.end();setIter++){
			VertexData *p = *setIter;
			vt += p->volume;
		}
		int j = 0;
		for(setIter = setPW.begin();setIter!=setPW.end();setIter++){
			vwell = new VertexWell;
			VertexData *p = *setIter;
			vwell->ID = p->ID;
			vwell->Qi = Qt*(p->volume/vt);
			vwell->vertex = p;
			pTMesh_data->productionWells_list[i].setValue(j++,0,vwell);
		}
//		cout << "Total volume for production well with flag " << mapIter->first << ": " << vt << endl;
//		cout << "#nodes " << setPW.size() << endl;
		i++;
	}

	i = 0;
	int numIW = mapIWells.size();											// total number of production wells
	pTMesh_data->setNumInjectionWells(numIW);
	pTMesh_data->injectionWells_list = new Matrix<VertexWell*>[numIW];			// allocate memory for a set of node well for each flag
	mapIter = mapIWells.begin();
	for(;mapIter!=mapIWells.end();mapIter++){
		double vt = .0;
		int flag = mapIter->first;
		pTMesh_data->setInjectionWellFlag(i,flag);
		well_info = MWells[flag];
		double Qt = well_info.flowRate;
		setWells setIW = mapIter->second;
		int num =  setIW.size();												// number of node for well flag
		pTMesh_data->injectionWells_list[i].allocateMemory(num,1);
		setWells::iterator setIter = setIW.begin();
		for(;setIter!=setIW.end();setIter++){
			VertexData *p = *setIter;
			vt += p->volume;
			//cout << "vt = " << vt << endl;
		}
		int j = 0;
		for(setIter = setIW.begin();setIter!=setIW.end();setIter++){
			vwell = new VertexWell;
			VertexData *p = *setIter;
			vwell->ID = p->ID;
			vwell->Qi = Qt*(p->volume/vt);
//			cout << "---------------------------------------------\n";
//			cout << "Node = " << p->ID << endl;
//			cout << "Qt = " << Qt << endl;
//			cout << "Qi = " << vwell->Qi << endl;
//			cout << "Vt = " << vt << endl;
//			cout << "Vi = " << p->volume << endl;
////			cout << "row = " << row << endl;
//			cout << "---------------------------------------------\n";
			vwell->vertex = p;
			pTMesh_data->injectionWells_list[i].setValue(j++,0,vwell);
		}
	//	STOP();
		cout << "Total volume for injection well with flag " << mapIter->first << ": " << vt << endl;
		cout << "#nodes " << setIW.size() << endl;
		i++;
	}

	int m,n;
//	double Qt = .0;
//	for (int i=0; i<numIW; i++){
//		pTMesh_data->injectionWells_list[i].getSize(m,n);
//		for (int j=0; j<m; j++){
//			VertexWell *vwell = pTMesh_data->injectionWells_list[i].getValue(j,0);
//			cout << "Injection well\tID " << vwell->ID << "\tQi = " << vwell->Qi << endl;
//			Qt += vwell->Qi;
//		}
//		cout << "total flow rate for well " << Qt << endl;
//	}
//	Qt = .0;
//	for (int i=0; i<numPW; i++){
//		pTMesh_data->productionWells_list[i].getSize(m,n);
//		for (int j=0; j<m; j++){
//			VertexWell *vwell = pTMesh_data->productionWells_list[i].getValue(j,0);
//			cout << "Production well\tID " << vwell->ID << "\tQi = " << vwell->Qi << endl;
//			Qt += vwell->Qi;
//		}
//		cout << "total flow rate for well " << Qt << endl;
//	}

	mapPWells.clear();
	mapIWells.clear();
	//throw Exception(__LINE__,__FILE__,"OK!\n");
}

bool SimulatorParameters::isInjectionWell(int flag) const{
	MWCIter mit = MWells.find( flag );
	return (mit!=MWells.end())?mit->second.isInjection:false;
}

bool SimulatorParameters::isProductionWell(int flag) const{
	MWCIter mit = MWells.find( flag );
	return (mit!=MWells.end())?mit->second.isProduction:false;
}

double SimulatorParameters::getFlowrateValue(int flag) const{
	MWCIter mit = MWells.find( flag );
	return mit->second.flowRate;
}

double SimulatorParameters::getTotalInjectionFlowRate() const{
	double Q = 0.0;
	for (MWCIter mwiter=MWells.begin(); mwiter!=MWells.end(); mwiter++){
		//	printf("well flag %d. Flow rate: %f.  Volume: %f\n",mwiter->first,mwiter->second.flowRate,mwiter->second.wellVolume);
		if ( isInjectionWell( mwiter->first ) ) Q += mwiter->second.flowRate;
	}
	return Q;
}

double SimulatorParameters::getWellVolume(int well_flag) const{
	MWCIter mit = MWells.find( well_flag );
	return mit->second.wellVolume;
}


bool SimulatorParameters::hasNodeWell(const int &flag){
	// if node has a prescribed saturation it's supposed that exist a well
	if ( mapSaturation.size() ){
		MapFlagIter MIter = mapSaturation.find(flag);
		return (MIter != mapSaturation.end());
	}
	else{
		isNodeFree(flag);
		return well;
	}
}
}
